import numpy as np
from scipy.sparse import csr_matrix


def mass_vector_matrix(vspace,mu):
    E = 2.0*mu
    mesh = vspace.mesh
    ldof = vspace.number_of_local_dofs()
    gdof = vspace.number_of_global_dofs()
    NC = mesh.number_of_cells()
    shape = (NC, ldof, ldof)
    cell2dof = vspace.cell_to_dof()
    I = np.broadcast_to(cell2dof[:, :, None], shape=shape)
    J = np.broadcast_to(cell2dof[:, None, :], shape=shape)

    bcs, ws = vspace.integrator.quadpts, vspace.integrator.weights

    rho = mesh.bc_to_point(bcs)[...,0] #(NQ,NC)
    phi = vspace.basis(bcs) #(NQ,1,ldof)
    

    M = np.einsum('i,ij,ijk,ijl,j->jkl',ws,1.0/rho,phi,E*phi,vspace.cellmeasure)


    M = csr_matrix((M.flat, (I.flat, J.flat)), shape=(gdof, gdof))

    return M

